The golden number seen in a mechanical oscillator

A seemingly ubiquitous irrational number often appearing in nature and in man-made things like structures, paintings, and physical systems, is the golden number. Here, we show that this astonishing number appears in the periodic solutions of an underactuated mass-spring oscillator driven by a nonlinear self-excitation. Specifically, by using the two-time scale perturbation method, it is analytically demonstrated that the golden number appears in the ratio of amplitudes, as well as in the oscillation frequency of the periodic solution, which is referred to as golden solution and, by applying the Poincaré method, it is demonstrated that this solution is asymptotically stable. Additionally, the analytic results are illustrated by means of numerical simulations and also, an experimental study is conducted.

(2) www.nature.com/scientificreports/ self-sustained oscillations in the system. Specifically, it is shown that the system has a stable periodic solution, which has the following properties: (1) the ratio of amplitudes corresponding to the actuated and non actuated oscillators is equal to the golden ratio, (2) the oscillation frequency of the periodic solution is equal to the natural frequency of the uncoupled oscillators scaled by the golden number. It is important to stress the fact that the aforementioned properties are independent of the intrinsic parameters of the oscillators. Furthermore, the occurrence of the golden number in the periodic solution is analytically demonstrated by using a perturbation method namely, the two-timescale method and also, analytic conditions for the stability of the periodic solution are derived by using the Poincaré method of perturbation. The obtained theoretical results are numerically illustrated and also, an experimental study is conducted in order to validate the analytic results. The experiments are performed using an electro-mechanical device. Ultimately, the obtained analytical, numerical, and experimental results presented here give clear evidence that the underactuated mechanical oscillator with weakly nonlinear self-excitation has "golden solutions": the golden number appears in the amplitude and frequency of the stable periodic solution.

System description and main result
Consider the mechanical system shown in Fig. 1, which consists of two identical masses of mass m (kg). The masses are coupled through a spring and also, one of the masses is attached to a fixed support via a spring. Both springs are identical and have stiffness coefficient k (N/m). Furthermore, only the first mass is actuated by the input signal U 1 . Consequently, the mechanical system is indeed an underactuated system because it has more degrees of freedom than actuators, see. e.g. Refs. 17,18 .
The dynamic behavior of the system is governed by the following set of equations where x 1 , x 2 ∈ R describe the positions of the masses. On the other hand, the following nonlinear term is applied to the actuated mass in order to induce a selfsustained oscillation where 0 < ǫ ≪ 1 is a small parameter and , H * ∈ R + are design parameters. Note that the nonlinear term (5) dissipates energy for ( 1 2 mẋ 2 1 + 1 2 kx 2 1 ) > H * , whereas it pumps energy into the system for ( 1 2 mẋ 2 1 + 1 2 kx 2 1 ) < H * , see e.g. Ref. 19 .
Determining the golden periodic solution in the system. In order to unveil the presence of the golden number in the amplitude and frequency of the periodic solution in system (4)-(5), we use the twotimescales method, also known as two-timing method, see e.g. Refs. 20-23 . First, let define the dimensionless time such that system (4)-(5) can be written in the form where the over-dots denote differentiation with respect to τ and (4) www.nature.com/scientificreports/ The next step is to look for solutions of (7) in the form Now, the two-timescale method is used in order to obtain explicit expressions for x 1 , x 2 , and a. In accordance with this well-known perturbation method and considering (9), the following ansatz is proposed where τ is the fast time and T = ǫτ is the slow time. Note that (10) is indeed and asymptotic approximation under the assumption that ǫ ≪ 1.
The first and second derivatives of (10) with respect to τ are given by where the primes indicate partial differentiation with respect to τ , that is x ′ = ∂x ∂τ . Replacing the approximated solution (10) and its derivatives (11) into the first equation of (7), and ignoring high order terms in ǫ (≥ ǫ 2 ) , yields to Before solving (12), it is convenient to re-scale the time according to such that Eq. (12) is rewritten as follows where now the primes indicate partial differentiation with respect to time τ.
By collecting terms with similar powers in ǫ , it is possible to obtain the following set of equations The general solution of (15) is where A(T) is the amplitude of the solution and θ(T) its phase. This result can be used for solving (16). First, note that where the over-dots denote differentiation with respect to T. Replacing (18) and (19) into (16) yields to with where P = p 1 ω (p 2 a 2 A 2 (ω 2 sin 2 (τ + θ) + cos 2 (τ + θ)) − H * )A sin (τ + θ). In order to continue with the analysis, it is worth to remember the following result from the theory of differential equations 23 : Let h(τ ) be a 2π-periodic function of τ . Let θ be any fixed real number. Then, the equation ẍ + x = h has a 2π-periodic solution x if and only if (9) x 1 (τ ) = −ax 2 (τ ), with a ∈ R + .
(17) x 01 (τ , T) = A(T) cos (τ + θ(T)), www.nature.com/scientificreports/ Since function h(τ ) given in (21) is indeed 2π-periodic, we can make use of Proposition 1. Then, replacing (21) in (22), and solving the integrals yields to From (23) it is possible to obtain the following differential equation which has the following equilibrium points The equilibrium point A ⋆ is unstable, whereas the equilibriums A ⋆ + || A ⋆ − are stables. Since A represents the amplitude of the periodic solution, see Eq. (17), we will focus, without lose of generality, on the positive equilibrium A ⋆ + . Then, by linearizing (25) around A ⋆ + it can be easily shown that the limit solution of (25) satisfies On the other hand, from (24) it follows that θ should satisfy which implies that θ is a constant, i.e.
with θ 0 ∈ R being an arbitrary constant. The next step is to replace (27) and (29) in the approximated solution (17), which results in Now, by substitution of (30) into the approximated solutions (10), and transforming from time τ back to τ using (13), yields to with A ⋆ + given in (27). Note that the approximated solutions (31)-(32) satisfy (9). The remaining task at this point is to find the value of a in (31)-(32). For this, the second derivative of (32) is computed, which results in Next, substitute (31)-(33) into the second equation of (7) to obtain where the argument of the cosine function is the same than in (32). In order to obtain the value of a, it is convenient to rewrite (34) in the form Clearly, for satisfying (35) for all time τ it is required that .
. www.nature.com/scientificreports/ Note that this polynomial in a is exactly the same than the polynomial from which the golden number is obtained, see Eq. (36). Therefore, as already mentioned in (3), the values of a satisfying (36) are However, since a should be a positive number, see (9), only the positive root of (36) is considered and therefore we have that a = ϕ:= 1+ √ 5 2 , i.e. the parameter a appearing in the amplitude and frequency of solutions (31)-(32) is the golden number.
Then, by substituting a = ϕ , (6), and (26) into (31)-(32), it follows that the approximated solution of system (4)-(5) is Finally, note that the above expressions can be further simplified by replacing p 1 and p 2 as given in (8) and ω given in (13). This yields Clearly, the golden number ϕ appears in the ratio of amplitudes and likewise, the oscillation frequency of the solution is equal to the natural frequency ω = √ k/m scaled by the golden number.

Stability of the golden solution.
In this part, the stability of the golden solutions (38)-(39) is investigated by using the Poincaré method of perturbation 24 . As a first step, system (7) is rewritten in the weakly nonlinear form It is important to note that in (42) we have added and subtracted the term ω 2 x 1 in the first equation of (7) and also the term ω 2 x 2 to the second equation. This may seem artificial, but in this way, we guarantee that for µ = 0 , the systems has the desired oscillation frequency. Furthermore, the obtained results will only be valid for that specific frequency.
Next, system (42) is diagonalized using the transformation where V is the matrix of eigenvectors associated to matrix A, see (43), and i = √ −1 . Then, in the new coordinates z ∈ R 4 , system (42) takes the form www.nature.com/scientificreports/ Note that for ǫ = 0 , the asymptotic solutions of (46) are where α i are arbitrary constants. However, for ǫ = 0 , the interest is in finding a unique set of α i 's such that the obtained solution exists and is asymptotically stable. This can be done by using Theorem 1 presented in Ref. 24 . That theorem establishes that the unique set of α i 's can be obtained by solving the following periodicity conditions where the functions P j , j = 1, . . . , 4 , are computed from where n 1 = n 3 = 1 , n 2 = n 4 = −1 , and f j (·) is the jth entry of vector V −1 F(Vz) , which is described in (46). Solving (51) yields where the lowercase p i 's are as given in (8).
After long but straightforward computations it can be shown that the values of α i 's satisfying (48)-(50) are where a is the solution of (36). Note that in the original coordinates x, the solutions (47), (56) exactly coincide with those already obtained in (38)-(39) with θ 0 = −π/2 . This can be easily shown by replacing (47), (56) into (45) with a = ϕ and using (13).
As a next step, the stability of solutions (47), (56) is investigated using a result based on the Poicaré method of perturbation. Specifically, the stability of the solutions is determined using again Theorem 1 in 24 . Following that theorem, yields to the characteristic polynomial where α * i , i = 1, 2, 3 , are the values of the α i 's given in (56). Then, according to Theorem 1 in Ref. 24 , the periodic solutions (47)-(56) are asymptotically stable if and only if all the roots of the characteristic polynomial (57) have negative real parts.
In order to determine this, let replace (48)-(50), with P j , j = 1, . . . , 4 , as given (52)-(55), into (57). This yields the polynomial where (46) www.nature.com/scientificreports/ At this point it is worth mentioning that a can only take two possible values, see Eq. (37). For any of these two values, the term (3a 2 − 3a + 2) in (60) satisfies (3a 2 − 3a + 2) = 5 , and likewise, the imaginary term in (61) vanishes. Consequently, the coefficients (60)-(61) are reduced to By straightforward computations, it is possible to show that all the roots of the characteristic polynomial (58) with coefficient a 2 given in (59), and coefficients a 1 and a 0 as given in (62) Summarizing, the above results provide analytic evidence that the golden number appears in the amplitude and frequency of the periodic solution of system (4) and that this solution is asymptotically stable. The analysis has been conducted using two perturbation methods namely, the two-timing and the Poincaré methods. The obtained results are formalized in the following theorem.
Theorem 1 Consider the dynamical system described by (4). Assume that the actuated part of the system is driven by the nonlinear term (5) with 0 < ǫ ≪ 1 . Then, the approximated periodic solutions ( O (ǫ) ) of system (4)- (5) have the form where ϕ ∈ R + is the golden number, ω = k m ∈ R + is the natural frequency, θ 0 ∈ R + is an arbitrary phase shift, and where H * is a design parameter of the nonlinear excitation (5) and k ∈ R + is the stiffness. Due to the fact that the golden number appears in the ratio of amplitudes and also in the oscillation frequency, solutions (65) are referred to as golden solutions. Furthermore, these solutions are asymptotically stable.

Remark 1
The golden number appears in the ratio of amplitudes of the approximated periodic solutions only if the oscillators are identical. Thus, the golden number can be used a measure of how identical or dissimilar the oscillators are.

The golden solutions and the energy in the system
This section highlights the importance of the golden solutions (65) from an energy viewpoint. For this, it is convenient to use the theory of passive systems and consider the following definition, which has been derived following 25 . Definition 1 Consider system (4) with input U 1 and output y =ẋ 1 , i.e. the output of the system is the velocity of the first oscillator. Then, the system is said to be lossless, i.e. with no energy dissipation, if there exists a continuously differentiable positive semidefinite function V(x) (called the storage function) such that For the present case, let consider the following storage function (60) a 1 = γ 2 π 2 a 4 (3a 2 − 3a + 2), 1 + a 2 )).
Remark 2 Note that there are two possible values of a for which (−a 2 + a + 1) = 0 in (70). However, as demonstrated in Theorem 1, the only stable periodic solution satisfying (9) is obtained for a = ϕ . Thus, the golden solution (65) is important from an energy viewpoint because for that particular solution there is no dissipation in the system.

Numerical study
In this section, the analytic results summarized in Theorem 1 are illustrated by means of numerical simulations. Hence, system (4)-(5) is numerically integrated using the Runge Kutta method implemented in Matlab with a step size of 1 × 10 −3 . The parameter values are summarized in Table 1 and the initial condition x 1 (0) is arbitrarily chosen as x 1 (0) = 5 × 10 −3 (m), and the remaining initial conditions are set to zero. Figure 2a shows the complete time series corresponding to the positions of the oscillators, denoted by x 1 (in blue) and x 2 (in green), whereas Figure 2b shows a zoom in the last two seconds of the simulation (steady behavior).
For the parameter values considered in this simulation, the expected values for the amplitudes, obtained using (65), (66), are x 1amp = ϕA * + = 0.0108 and x 2amp = A * + = 0.0067 . These values are indicated by the horizontal blue and green dashed lines respectively, in Fig. 2b. On the other hand, the amplitudes obtained from the simulation (solid blue and green signals) are x 1amp = ϕA * + = 0.01081 and x 2amp = A * + = 0.00668 . Clearly, there is a good agreement between the theoretical and numerical results. Consequently, these numerical results confirm that the ratio of amplitudes coincides indeed with the golden number. This can be further verified from Fig. 3, which shows the projections of the solutions of system (4)-(5) onto the ( x 2 , x 1 )-plane, see the blue line. The dotted red line is just a reference line with slope −ϕ.
Finally, Fig. 4 shows the single-sided amplitude spectrum for solutions x 1 and x 2 . Clearly, the oscillation frequency of the solutions is at ϕω = ϕ k m = 20.9356 (rad). Note that the oscillation frequency is indeed equal to the natural frequency scaled by the golden number, just as indicated by Theorem 1.
In summary, the numerical results presented here give clear evidence that the golden number appears in the amplitude and in the frequency of the periodic solutions of system (4)- (5).

Experiments
This section presents and experimental study in order to validate the theoretical and numerical results obtained in previous sections.
The experiments are conducted in the rectilinear plant ECP model 210, shown in Fig. 5, which consists of two mass-spring oscillators. The positions of the masses are obtained by means of encoders and one of the masses can be directly actuated by means of a servo actuator. The parameter values for the system and for the nonlinear term (5) are as given in Table 1 and a data acquisition card is used for measuring the signals from the encoders and also, for sending the nonlinear term (5) to the servo actuator. Furthermore, the velocity of mass one, which is required in (5), is obtained via a Luenberger observer.
The time series obtained from the experiment for x 1 (blue) and x 2 (green), are shown in Fig. 6. The amplitude of x 1 is x 1amp = 0.0105 (m), whereas the amplitude for x 2 is x 2amp = 0.0065 (m).
Note that the ratio of amplitudes x 1amp /x 2amp = 1.6180 is very close to the golden number ϕ = 1.6180 . . . , there is just a tiny error of 0.8%.
(70) V = (−a 2 + a + 1)kx 2ẋ2 + U 1 y. Table 1. Parameter values of system (4)- (5). These parameter values have been obtained partly from the experimental setup ECP MODEL 210 discussed in the next section. www.nature.com/scientificreports/ To further illustrate this, the measured time series are projected onto the ( x 2 , x 1 )-plane. This is shown in Fig. 7, where the blue line is the experimental result and the red line is an auxiliary line with slope equal to −ϕ . Note however, that the blue curve is not a perfect straight line due to the nonlinear effects like dry friction and the unavoidable differences between the oscillators. Nevertheless, the observed agreement is remarkable. Now, for obtaining the frequency from the experimental data, it is convenient to take into account the following. Solutions x i , i = 1, 2 are very well approximated by x i (t) = A i sinωt , with ω the unknown frequency to be determined. Furthermore, the derivative of  www.nature.com/scientificreports/ A i is already known, see Fig. 6, then the unknown frequency can be obtained by dividing the amplitude of the derivative between the amplitude of the signal, i.e. ω = B i /A i . In order to get B i , the approximated derivative D x i of x i , also known as dirty derivative 26 , is computed as follows for n = 1, . . . , m − 1 , where m is the length of the signal. Note that this approximation holds for t n+1 − t n ≪ 1 . Then, the value of B i is obtained from B i = max(D x i ) . Using x 1 for computing the approximated derivative, it follows that B 1 = 0.2166 . Consequently, the approximated oscillation frequency is ω = 20.632 (rad). On the other hand, the expected oscillation frequency, see Theorem 1, is ϕω = ϕ k m = 20.936 (rad). Then there is an error of 1.45% between the experimental and theoretical results, which is acceptable given the fact that the theoretical analysis has been conducted under the assumption of identical oscillators, which in practice is impossible to achieve.  www.nature.com/scientificreports/ Additionally, in order to illustrate Remark 1, we have repeated the above experiment but this time the oscillators have different masses. Specifically, an extra mass is added to oscillator 1. The obtained results are shown in Fig. 8. The blue curve corresponds to the case that an extra mass of 0.250 (kg) is added, whereas the orange curve corresponds to an extra mass of 0.750 (kg). From this figure it is clear to see that as the difference in mass increases, the ratio of amplitudes deviates more and more from the golden number (yellow line).
Ultimately, the experimental results presented here give strong support to the theoretical and numerical results presented in previous sections.

Discussion and conclusions
We have provided analytical, numerical, and experimental evidence that the golden number appears in the amplitude and frequency of the asymptotically stable periodic solution corresponding to an underactuated mechanical system with weakly nonlinear self-excitation.
For the case of identical oscillators, the absolute ratio of the amplitudes coincides exactly with the golden number and likewise, the oscillation frequency is given by the natural frequency of the oscillators scaled by the golden number. Interestingly, this is true independently of the value of the intrinsic parameters of the system like mass and stiffness, as long as the assumption of weak nonlinearity ( ǫ ≪ 1 ) is satisfied. Furthermore, it has been shown that if the ratio of amplitudes of the solutions of the coupled oscillators coincide with the golden number, then the system is lossless, i.e. there is no energy dissipation in the system.
On the other hand, the experimental results presented here have shown a remarkable approximation to the theoretical results besides the fact that the latter have been obtained assuming identical oscillators, whereas  www.nature.com/scientificreports/ in the former this assumption is obviously never satisfied. In fact, while performing the experiments, we have noticed that for the system under consideration, the golden number can be used as a measure of how identical or dissimilar the oscillators are. Specifically, if the absolute value of the ratio of amplitudes is very close to the golden number then it is possible to conclude that the oscillators are practically identical without doing any extra computation. Ultimately, based on the results reported in this manuscript, we claim that the mechanical system discussed here can be added to the list of systems/examples where the golden number emerges. Furthermore, since there exists a well-known relationship between mechanical and electrical circuits, see e.g. Ref. 27 , we believe that the results presented here may also be applicable to certain electronics circuits, like for example capacitively coupled LC circuits.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.